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We propose an efficient and fast numerical algorithm of finding a stationary solution 
of large systems of aggregation-fragmentation equations of Smoluchowski type for con¬ 
centrations of reacting particles. This method is applicable when the stationary concen¬ 
trations steeply decreases with increasing aggregate size, which is fulfilled for the most 
important cases. We show that under rather mild restrictions, imposed on the kernel of 
the Smoluchowki equation, the following numerical procedure may be used: First, a com¬ 
plete solution for a relatively small number of equations (a ’’seed system”) is generated 
and then the result is exploited in a fast iterative scheme. In this way the new algorithm 
allows to obtain a steady-state solution for rather large systems of equations, by orders 
of magnitude faster than the standard schemes. 

Keywords : kinetic equations; numerical algorithm; aggregation 


1. Introduction 

Aggregation and fragmentation are ubiquitous in nature and many industrial pro¬ 
cesses are based on these physical phenomena. As a representative example one can 
mention reversible polymerization in solutions, coagulation in colloidal solutions, 
formation of proto-planets in dust clouds in the inter-galaxy space, etc. The size 
distribution of the particles in planetary rings, like Saturn rings, is also determined 
by a subtle balance between aggregation and fragmentation of the rings objects. 
Aggregation processes play an important role in living systems: Coagulation of pri¬ 
ons takes place in numerous prion-related diseases; moreover, the blood clotting 

pm 

may be mentioned as the most vivid example from the everyday life 1 |0 1 For 
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space homogeneous systems, addressed in here the above processes are described 
by Smolucliowski-like equations. These quantify the change of the size of particles 
due to formation of joint aggregates from colliding entities, or due to fragmenta¬ 
tion of the aggregates. While the former process is necessarily binary, the latter 
one may be either binary, when particles break upon a collision or unary, when a 
decay of an aggregate into smaller pieces takes place spontaneously. Binary frag¬ 
mentation occurs due to inter-particle impacts, when part of the kinetic energy of 
colliding objects (for instance, of ice grains in planetary rings) exceeds some frag¬ 
mentation threshold. Unary decomposition may occur due to thermal fluctuations, 
for example, in the case of decomposition of polymers®, or in Saturn rings due to 
meteoroid bombardment. In both cases of spontaneous or collisional fragmentation 
the distribution of debris has some particular form, depending on the nature of the 
fragmentation process. 

In the present study only binary collisions are addressed, that is, we do not con¬ 
sider inter-particle reactions which involve three or more particles. This is justified 
for the most important applications; it is a rather seldom process, when an aggre¬ 
gation or fragmentation process may occur if more than two particles are involved. 

Below we will study the following fragmentation processes: (i) spontaneous de¬ 
composition due to thermal fluctuations (unary process) (ii) fragmentation due to 
inter-particle impact, when part of the kinetic energy of colliding objects (for in¬ 
stance ice grains in planetary rings) transforms into the surface energy of debris. 
While in the former case fragmentation occurs as a random event, the latter one is 
possible if the kinetic energy of particles exceeds some threshold 

We assume that all aggregates in the system are comprised of particles of minimal 
size, monomers, which do not suffer further splitting 1 A The physical nature of 
the monomers may be very different, depending on the nature of the system. For 
instance, in polymer solution the role of monomers play some functional chemical 
groups, which can associate into larger molecules. In planetary rings, the role of 
monomers play ice grains of a minimal size observed in these system; such grains, 
however, may be formed from still smaller particles. 

Let rife (t) be concentration (a number of particles in a unit volume) of aggregates 
comprised of k monomers, of mass mi and radius ri, at time t. The Smoluchowski 
equations describe time evolution of the concentrations rife(t). As a system evolves, 
larger and larger aggregates emerge, hence, strictly speaking the system of Smolu¬ 
chowski equations must be infinite. Obviously, a direct application of standard nu¬ 
merical schemes to solve such infinite systems is not possible, so we need to develop 
a more advanced approach. In practice, only a finite number of equations may be 
numerically solved, therefore a problem arises, how to mimic a behavior of an infi¬ 
nite system with a use of a finite subset; moreover, the basic physical principles, like 
conservation of mass should not be violated. Next problem is the large number of 
equations, needed for an adequate modeling, even if a finite system is used. Indeed, 
the number of different species (that is, aggregates of different sizes) may be very 
large in some systems; it is equal, for instance, to 10 9 — 10 12 for planetary rings, 
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of the possible solutions is the application of the coarse-graining, where concentra¬ 
tions of aggregates with close number of monomers are grouped together; the larger 
is the size of an aggregate, the larger is the respective group®. Unfortunately, this 
approach is not always possible, especially in the case when a sharp cross-over from 
one functional form of size distribution to another one takes place. As we will show 
below such cross-over is indeed observed for some of the studied systems. The re are 
other methods to handle large systems of Smoluchowski equations, e.g. 
these methods, are focused, however, on the time-dependent solu tions and do not 
take precautions against violation of the mass conservation 10 H. The latter is of 
a primary importance for systems considered here. In the present study we analyze 
the stationary solutions of large systems of Smoluchowski equations and develop 
an efficient numerical algorithm. It provides a fast solution for the particle size dis¬ 
tribution, even if the concentrations of different species differ by 10-12 orders of 
magnitude. The rest of the paper is organized as follows. In the next Sec. II we 
consider the aggregation-fragmentation equations of the Smoluchowski’s type. In 
the Sec. Ill we discuss the most widely used models for the rate coefficients of the 
Smoluchowski equation and provide the existent exact solution for some particular 
model. In Sec. IV we discuss in detail the new numerical approach and present 
simulation results obtained with the new method. Here we compare the efficiency of 
our method with the traditional one. In the last Sec. V we summarize our findings. 

2. Aggregation-fragmentation models 

To introduce the notations we start with the standard Smoluchowski aggregation 
equations. The aggregation process may be symbolically written as 



In what follows we will consider space uniform systems, that is, the local concen¬ 
trations rii(r,t) are assumed to be uniform trough the system, = rii(t). We 

also introduce the kinetic coefficients Cij, which give the production rate of aggre¬ 
gates of size i + j from particles of size i and j. C\j quantifies the number of such 
aggregates that appear in a unit volume during a unit time. In these notations the 
kinetic equations read: 



( 1 ) 


i+j=k 


The equations m with k = 1,2,... are the standard Smoluchowski equations. The 
first term in the r.h.s. of the above equation describes the rate at which aggregates 
of size k are formed from particles of size i and j. The summation extends over all 


a Indeed, while the pr imary particles have the radius of ~ 1 cm, the maximal aggregate size could 
be about 5-10 meters mm; 
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i and j with i + j = k, and the factor ^ prevents double counting. The second term 
gives the rate at which the particles of size k disappear due to aggregation with 
other particles of any size i. 

Fragmentation is the opposite process when an aggregate splits into smaller 
pieces. It may happen spontaneously, due to thermal fluctuations, when, for instance 
a particle of size k breaks into / fragments of size * 1 , 12 , ■ ■ - h, 

[*] —* [*1] + M + • ■ • [n] 

where *i + ii + ... ii = k, due to the mass conservation. 

The fragmentation is also possible due to collisions of aggregates. In this case the 
kinetic energy of their relative motion transforms into kinetic energy and surface 
energy of debris. This process may be symbolically written as 

\ k ] + bl —* [* 1 ] + N + ■ ■ ■ [ii\ 

where again *i + *2 + • ■ • ii = k + j due to the mass conservation. Let us consider a 
simplified collision model, assuming that particles break completely in a disruptive 
collision into monomers of mass mi. The kinetic equations, describing all the above 
processes of aggregation as well as spontaneous and collision fragmentation read, 

^ = 2 ^ . CijTliTlj 'y ^ (Cik T Aik) BiTlk AkTlk (2) 

i+j=k i> 1 

-^- = ~ ni H Cl i n i + ni j A V n j + \ Yh Ai i (* + fi nin i + H j A i n d ■ 

3> 1 j> 2 i,j> 2 j> 2 

Here Ak is the rate of spontaneous fragmentation and Aij - the rate of collisional 
fragmentation. Note, that similar equations have been studied in the context of of 
rain drop formation . 

The dependence of aggregation and fragmentation rate coefficients on the parti¬ 
cles size is determined by the nature of the process. For example, the rate coefficients 
for the diffusion-limited aggregation read ^ 

Cij = Co(i 1/3 + j 1/3 )(i~ 1/3 + j~ 1/3 ) (3) 

where the constant Co depends on the characteristic size of agglomerates and their 
char acteri stic mobility. In the case of ballistic aggregation such dependence has the 

f.in !4JJ1 

c„ = c;(.‘''“+j‘«) 2 (r‘+r I ) ,/2 , (4) 

where again, C' 0 is determined by the characteristic size and velocity of the particles. 
Noteworthy, that the rate coefficients |3]) and © are homogeneous functions of 
masses of colliding particles nii ~ i and rrij ~ j. Sometimes simplified rate kernels 
are considered 


Ch = Co (ijY . 


( 5 ) 
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It may be shown that for some important applications the collision frag mentation 


coefficients A i: j are proportional to the aggregation coefficients, e.g. 




A jj — \Cij , 


( 6 ) 


where A is a constant that characterizes the ratio of the disruptive and aggregative 
collisions. Usually, it is rather small, A <C 1. For the case of spontaneous fragmen¬ 
tation the following model may be applied, 


At. = vk H 


( 7 ) 


where the exponent 6 depends on the dimension of aggregates (which may be, 
generally fractal) and the average coordination number of the primary particles 
packing. 

In what follows we describe the fast numerical algorithm of finding the steady- 
state solution of the large sets of Smoluchowski-type equations ©. It is also impor¬ 
tant to check the accuracy of the novel approach. Therefore we give here, without 
derivation, the exact steady-state solution to these equations for the particular case 
of constant rate coefficients, Cij = Co, Aij = ACo with A = const and Aj~ = v , that 
is, for fjt = 6 = (PUH 


n k = 


n\ = 


M 0 = 


2ni 


i r(fc-|) 

yj 47r r (fc + 1) V, (1 + A)M 0 + v 
M (AM 0 + v) 
v+{l + A )M 0 

A M-v+{{y- AM) 2 + 2 vM (1 + 2A)) 
1 + 2A 


[(1 + A) M 0 + v\ 


1/2 


( 8 ) 


Here F(a;) is the gamma function, M = ^ n k is the total mass and M 0 = 

Sfci n k ~ the total number of the aggregates; without the loss of generality we 
take Co = Ao = 1. The derivation details for the above solution will be given in the 
forthcoming publication ^ ^. 


3. Fast numerical algorithm 
3.1. Description of the algorithm 

Solving the system of Smolucliovski equation one encounters a number of technical 
problems. The first problem in the numerical analysis of the system of infinite 
number of rate equations © is the conservation of mass of the constituents. Indeed, 
in any real simulation one can handle only a finite number of equations say TV, 
which describe evolution of particles of size 1, 2, ...TV (a particle of size k has 
mass TO*, = m\k). These equations have both aggregation and fragmentation terms. 
In particular, they have a term which describes aggregation of particles of size 
i < N and j < TV, resulting in an aggregate of size % + j > TV. Since the system 
of TV equations does not account for particles larger that TV, such processes would 
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effectively lead to the leak of particles’ mass and violation of the mass conservation. 
To preserve the mass conservation we consider the model, where all collisions with 
i + j > N are disruptive, that is, such collisions produce i + j monomers. We have 
checked that this additional assumption does not lead to any noticeable distortion 
of the numerical solution of the rate equations n k , provided N is large enough and 
k is smaller than some fraction of N. 

The most significant problem is to efficiently handle a large number of equa¬ 
tions, say up to ~ 10 10 . In principle, this may be done by a coarse-graining, that 
is, by grouping all concentration n k ~ n k +i into coarse-grained variables fix with 
increasing l as k grows In the case of interest, however, we have a drastic vari¬ 
ation of the functional dependence of n k {k ), which changes from a power-law to 
the exponential decay. This hinders an effective application of the coarse-graining 
and we need to keep explicitly all individual concentrations. Hence we have to work 
with a large number of equations, which is computationally costly. To speed up the 
computations we have developed the following recursive procedure for the solution 
of system of equations <[2|). Taking into account, that we search for a stationary 
solution, drik+i/dt = 0, we obtain for the number density n k +i'- 


1 

2(1 +A) 


E CijUiUj 

i-\-j=k -\-1 


N 

^ ^ Ci k-\-l^tittk-\-l 
i =1 


Ai. 


1 + A 


n k = 0. 


(9) 


The first sum in Eq. ([9]) contains only rii with i < fc, while we write the second sum 
as 


N k N 

^ ' Ci k+l’kli’klk+l — ^k- f-1 y ) Ci k -\-\Tli + Ck+l k+l^k+l ~b ^k+ 1 y ' Ci k+l^li (.10) 
2— 1 2=1 2=fc+l 

Now we use the properties of the kinetic kernel C^ and the steady-state distribution 
Tik = n-k{mk), which we assume to be a decreasing function of k. Namely, we assume 
that the coefficients Cij remain constant or increase with i and j at a smaller rate 
than the rate at which n k decreases with k. That is, we assume that for k 3> 1 the 
following condition holds true: 

k N 

y ] Ci k -\-1 klj +> y ( C i k-\-lTli • (H) 

2=1 2=fc+l 

This allows to neglect the last sum in Eq. (HI and obtain the quadratic equation 
for n k+ 1 : 


Ck+l k+lKk+l + n k+l 


A 


■fc+i 


1 + A 




'2 fc+1^-2 


2=1 


E 

i-\-j=k -\-1 


CijTliTlj 

2(1 +A) 


= 0. 


Solving the above equation and choosing the positive root, we arrive at the recurrent 
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relation for the concentrations 



( 12 ) 


2 



k 


Fig. 1. Equilibrium size distribution of particles for Cij = 1, Aij = A, Ak = v- The total 
number of equation is N = 16000 and the size of the ’’seed system” is K = 1500 (see the text 
for the explanations). A good agreement between the numerical (solid lines) and exact solution, 
Eq. <0 (dashed lines) is observed. A slight deviation from the exact solution for large k may be 
attributed to the computational errors, when very small numbers are handled. 

Using the recurrent relation m one can significantly accelerate computations. 
This may be done as follows: First, we solve explicitly the system of K rate equa¬ 
tions (0 for K <C N, choosing the value of K to fulfill the condition (flTl) . Next, 
the concentrations rii with K < i < N may be straightforwardly obtained from 
the recurrent relation (1121) . In other words, our new algorithm uses the following 
trick: First we ’’exactly” solve a relatively small ’’seed system” and than apply the 
recurrence to obtain the solution for the whole system. 

Solving numerically the rate equations with different kernels directly, and with 
the use of our recurrent method, we proved the efficiency and accuracy of the 
above accelerating approach. The computation accuracy may be checked by a direct 
comparison of the numerical and exact solutions. The results of numerical solution 
of system of equations 0 is given in Fig. |T] A good agreement with the theory is 
observed. 
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3.2. Analysis of the algorithm efficiency 

The efficiency of our new algorithm can be evaluated in the following way. The 
direct solution of the system of equations requires ~ TV 2 operations at each time 
step. Let T be the number of time steps, needed to arrive at the steady state. 
Since the relaxation time depends only on the physical parameters of the system, T 
should not be sensitive to the number of the equations TV. Hence, it is reasonable to 
approximate it by a constant, T(TV) = const T?> 1. Then the number of operations 
t(TV), required to find a steady-state solution of a system of TV equations reads, 
t(N) = TN 2 . The number of operations t(N) refers to the standard computation 
scheme. 

Consider now the number of operations needed in our new approach for a system 
of TV equation. Let we start with the ’’seed system” of K equations. To find a 
steady-state solution for the first AT concentrations we need TK 2 operations. Then 
to obtain the steady state solution of the (K + l)-st concentration nx+i, using the 
fist I\ concentrations n±, n 2 , ■ ■ ■, n k, one needs K + 1 operation, see Eq. m- For 
the next (K + 2)-th concentration, (K + 2) operations are required, etc. Hence, to 
find the rest (TV — AT) steady-state concentrations one has to perform 

(AT + 1) + (K + 2)... + TV = i(TV - A")(TV + AT - 1) ~ i(TV 2 - k 2 ) 

operations, so that TAT 2 + (TV 2 — k 2 )/ 2 operations in total are needed to obtain the 
solution within the new approach. Then the ratio of computation times, required 
for the direct (standard) and new algorithm reads: 

t iec _ TK 2 + (N 2 — K 2 )/2 _ A' 2 1 AT 2 K 2 

f stan ~ TN 2 TV 2 + 2 T 2N 2 T ~ TV 2 ’ ^ ] 

where we use the fact that TV K and that the number of time steps, T T§> 1, 
which refers to the relaxation time, is always in practice much larger that the ratio 
N 2 /K 2 . 


Table 1. Comparison of the CPU time of the direct (standard) 
and new algorithms for the solution of the Smoluchowski equations 
for the case of collisional fragmentation for different number of 
equations. Cij = 1, Aij = A and A^ = 0 


A 

K 

N 

tree (a.e.) 

^stan (a.e.) 

^stan/^rec 

0.1 

10 3 

5 • 10 3 

37 

1787 

48 

0.05 

To 3 

To 3 

38 

10800 

284 

0.01 

1.5 ■ 10 3 

io 3 

80 

10090 

126 


In Table [T| we present the comparison of the efficiency (in terms of the CPU 
time) of our new algorithm and of the standard one; as it may be seen from the 
table, the new algorithm is indeed much faster and efficient. Note that the significant 
acceleration in the systems solutions is obtained without any loss in the computation 
accuracy. 
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4. Conclusion 

We perform a numerical analysis of Smoluchowski-like aggregation-fragmentation 
equations to find a steady-state size distribution of particles. We address very large 
systems of equations as required in some important applications. We develop a new 
numerical algorithm which yields a significant acceleration of simulations without 
the loss of its accuracy. We compare our method to the standard one and confirm 
its impressive superiority in speed. The developed algorithm may be exploited not 
only for the Smoluchowski-like aggregation-fragmentation equations but for other 
set of equations which have the structure similar to that addressed in our study. 
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